Vortex line in spin-orbit coupled atomic Fermi gases 
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It has recently been shown that the spin-orbit coupling gives rise to topologically-nontrivial and 
thermodynamically-stable gapless superfluid phases when the pseudo-spin populations of an atomic Fermi gas 
is imbalanced, with the possibility of featuring Majorana zero-energy quasiparticles. In this paper, we consider 
a Rashba-type spin-orbit coupling, and use the Bogoliubov-de Gennes formalism to analyze a single vortex line 
along a finite cylinder with a periodic boundary condition. We show that the signatures for the appearance of 
core- and edge-bound states can be directly found in the density of single-particle states and particle-current 
density. In particular, we find that the pseudo-spin components counterflow near the edge of the cylinder, the 
strength of which increases with increasing spin-orbit coupling. 

PACS numbers: 05.3Q.Fk, 03.75.Ss, 03.75.Hh 
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I. INTRODUCTION 

Following the recent experimental success with artificial 
gauge fields and spin-orbit coupled atomic Bose gases fHO], 
there has been increasing theoretical interest in studying spin- 
orbit coupled atomic Fermi gases with balanced or imbal- 
anced populations, at zero or finite temperatures, in two or 
three dimensions, etc. iH-[T^]. The main motivation behind 
these works is that the spin-orbit coupled atomic Fermi gases 
are ideal systems for studying topologically-nontrivial super- 
fluid phases |H,0-01 5 with the possibility of featuring Majo- 
rana zero-energy bound states for which the associated quasi- 
particle operators are self-Hermitian. This means that a zero- 
energy Majorana quasiparticle is its own anti-quasiparticle. 
Although these quasiparticles are predicted to appear in low- 
dimensional strongly-correlated systems in various fields of 
physics, including the fractional quantum Hall systems 11711 . 
chiral two-dimensional superconductors lfl8l[l9ll . chiral two- 
dimensional p-wave superfluids ll20l |2lll . three-dimensional 
topological insulator-superconductor heterostructures ll22tl 
one-dimensional nanowires J23l l24ll . spin-orbit coupled 
semiconductor-superconductor heterostructures l25ll26ll . etc., 
it has proved to be very difficult to realize them in these sys- 
tems. Given that the cold atom systems offer unprecedented 
control in comparison to condensed matter ones, there is a 
good chance of creating and observing Majorana quasiparti- 
cles with atomic systems in the near future. 

The first step in searching for the Majorana quasiparticles 
with spin-orbit coupled Fermi gases is to understand the phase 
diagram of these systems, which has recently been worked 
out within the mean-field approximation d, For in- 

stance, the ground-state phase diagram of a Rashba-type spin- 
orbit coupling at unitarity, i.e. when the two-body scatter- 
ing length a s between pseudo-spin components in vacuum di- 
verges, is illustrated in Fig. Q] There are three phases in the 
phase diagram While the normal (N) phase is char- 

acterized by a vanishing superfluid order parameter, the uni- 
form superfluid and nonuniform superfluid, e.g. phase separa- 
tion (PS), are distinguished by their thermodynamic stability 
when the order parameter is nonzero. Furthermore, in addi- 
tion to the topologically-trivial gapped superfluid (SF) phase, 
the gapless superfluid (GSF) phase can be distinguished by 



the momentum-space topology of its excitations. Depending 
on the number of zero-quasiparticle excitation energy regions 
in momentum space, there are two topologically-distinct gap- 
less phases. For the Rashba-type spin-orbit coupling shown 
here, while GSF(II) has four zero-energy points, GSF(I) has 
only two. 

The phase diagram illustrates that the spin-orbit coupling 
counteracts the population imbalance, and that this com- 
petition tends to stabilize the GSF phase against PS. The 
anisotropic nature of the spin-orbit coupling (in momentum 
space) is also found to stabilize exotic superfluid phases. For 
instance, in sharp contrast to the a — case where only the 
gapless superfluid phase supports population imbalance, both 
the gapless and gapped superfluid phases are found to support 
population imbalance. Although Rashba-type spin-orbit cou- 
pling is considered in Fig.Q] the topological structure shown 
here is quite robust against the effects of anisotropic spin-orbit 
couplings 0]. 
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FIG. 1: (Color online) The mean-field ground-state phase diagram 
of a Rashba-type spin-orbit coupled Fermi gas is shown as a function 
of the population-imbalance parameter P = (Nf — N±)/N and the 
strength of the spin-orbit coupling a at unitarity, i.e. when \a s \ — s- 
00. Here, N, SF, GSF, and PS correspond to normal, topologically- 
trivial gapped superfluid, topologically-nontrivial gapless superfluid, 
and phase separation, respectively. (Adapted from Ref. @].) 

Since Majorana quasiparticles appear in the presence of 
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topological defects, e.g. vortices, domain walls, boundaries 
between bulk phases, etc., one of the exciting experimental 
directions with spin-orbit coupled Fermi gases may be to cre- 
ate and observe a Majorana zero-energy quasiparticle bound 
to a vortex. For this purpose, in this paper, we consider a 
Rashba-type spin-orbit coupling, and use the Bogoliubov-de 
Gennes (BdG) formalism to analyze a single vortex line along 
a finite cylinder with a periodic boundary condition. We find 
signatures for the appearance of core- and edge-bound states 
in various observables, most notable of which is the counter- 
flow of pseudo-spin particle-current densities near the edge of 
the cylinder. The rest of the paper is organized as follows. 
First, we generalize the BdG formalism to spin-orbit coupled 
Fermi gases in Sec. HJ and then derive the self-consistency 
equations for a single vortex line in Sec. [TTTJ We present the 
numerical solutions in Sec. lIVl where we discuss the effects of 
spin-orbit coupling on the superfiuid order parameter, particle 
density, energy spectrum, bound-state wave functions, single- 
particle density of states, and probability-current density. Our 
conclusions are briefly summarized in Sec. [V] 



II. BOGOLIUBOV-DE GENNES FORMALISM 

The results mentioned above are obtained by solving the 
self-consistent BdG equations, suitably generalized to spin- 
orbit coupled Fermi gases. Before presenting our numerical 
results, let us first present the theoretical framework of this 
formalism. 



A. Hamiltonian 

In this paper, we use the mean-field Hamiltonian density (in 
units of H= k B = 1), H(r) = J2 a <T , V>£(r)i^(r)Vv(r) + 

A(r)'0|(r)'0|(r) + A*(r)^(r)^-|-(r), describing two- 
component Fermi gases with attractive and short-range 
interactions. Here, the operators V4-( r ) ar, d create and 

annihilate a pseudo-spin a fermion at position r, respectively, 
and A(r) is the mean-field superfiuid order parameter. 
Furthermore, the operator K aa (r) = — V 2 /(2M CT ) — 
is the kinetic energy, where M a is the mass and [i a is 
the chemical potential of a fermions, and the operator 
K^i(r) = ifL(r) = a(p y + ip x ) is the spin-orbit coupling, 
where a > is its strength and pj = —id/dj is the momen- 
tum operator. Although we keep the formalism quite general, 
we present numerical calculations only for the mass-balanced 
Fermi gases with = A/j, = M. 

In the presence of a spin-orbit coupling, this Hamilto- 
nian can be diagonalized via the generalized Bogoliubov- 
Valatin transformations, Vv( r ) = J2 n a 1 ( r )7«^' + 
v° a ' (r)7^ l(T /], where u°°' (r) and u£ CT '(r) are the ampli- 
tudes, and 7^ a and 7„. (T are the operators corresponding to 
the creation and annihilation of pseudo-spin a quasiparticles, 
respectively. The resultant BdG equation can be written as 



H(r)rt (r) 



H(r) = 



( r ). where 



K n (r) 


A*(r) 



-A*(r) 






-A(r) 
-K* n (v) 
-Kl t (r) 



A(r) 


-jq^(r) 



(1) 



is the Hamiltonian matrix given in the ip° (r) = 
[ul a \r),ui°\r),vl a '(r),v} l a '(r)] T basis, and e£ > are 
the energy eigenvalues. Since the BdG equations do not de- 
pend on a', it is sufficient to solve only for it^(r) = u™ (r), 
<(r)=tC'(r)and e n = fg. 



B. Self-consistency equations 

Using the Bogoliubov-Valatin transformations, the mean- 
field superfiuid order parameter A(r) = g(tpf(r)ipi(r)), 
where g > is the strength of the attractive interac- 
tion between f and I fermions, and (■ ■ ■ ) is the ther- 
mal average, becomes A(r) = gJ2 n i u n( r ) v i*( r )f( £ n) + 
4(r)4* (r)/(-£„). Here, f(x) = l/(e^ T + 1) is the 
Fermi function and T is the temperature. As usual, we relate 
the interaction strength g to the two-body scattering length 
a s between an f and a ! fermion in vacuum via the rela- 
tion, 1/g = -M r V/(4na s ) + ^ k l/(£ k ,t + £ia), where 
M r = 2M^-Mi/(Mf + M|J is twice the reduced mass of 
an j" and a J, fermion, V is the volume of the sample and 
ffk.er = k 2 /(2M a ) is the kinetic energy. This leads to g = 
4ir 2 a s /[2M r a s iJ2M r £ c — M r ir], where e c is the energy cut- 
off used in the k-space integration (to be specified below in 
SecUVj. 

To determine fi a , the order parameter equation has to be 
solved self-consistently with the number equations N a = 
J rfrji„(r), where ^^(r) = (?/4.(r)?/v(r)) is the local density 
of a fermions. Using the Bogoliubov-Valatin transformations, 
we obtain n CT (r) - EJKMI'/M + |<(r)| 2 f(-e n )}. 
Having generalized the BdG formalism to spin-orbit coupled 
Fermi gases, next we apply it for a single vortex line. 



III. SINGLE VORTEX LINE 

In particular, we consider a single vortex line positioned 
along a finite cylinder of radius R and length L, and with 
a periodic boundary condition in the z direction, in such 
a way that the order parameter can be written as A(r) = 
A(r)e~ 10 , where r and are the cylindrical coordinates r = 
(r, 9, z) flln . Note in this coordinate system that the spin-orbit 
coupling term becomes A"-j-j,(r) = e~ l0 [d/dr — (i/r)d/d8], 
showing that the single vortex line has rotational invariance 
around the z axis, so that the solutions of the BdG equation 
have a well-defined planar angular momentum m, i.e. m is 
a good quantum number. In addition, the system is assumed 
to have translational invariance along the z direction, i.e. p z 
momentum is also a good quantum number. 
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Thus, for a singly-quantized vortex line considered in 
this paper, we may choose the normalized wave func- 
tions as u£(r) = iij m5 (r ) e lm9 e lks z j\J 2nL and uj(r) = 
v} lms {r)e li - m+2 ' ,e e lksZ /^/2-kL for the t components, and 
ui(r) = ui ms (r)e l{m+1 ^e ik ^/V2^L and u£(r) = 
vi ms {r)e t{m+r)e e tk ° z /V2^L for the | ones. Here, k s = 
2irs/L is the wave vector along the z direction with 
s = 0, ±1,±2, .... This particular choice (which is not 
unique J25l HfJ) allow us to decouple the BdG equations 
into independent subspaces of (m, s) sectors. We further 
project the radial wave functions onto a set of Bessel func- 
tions normalized in a disc of radius R [27], i.e. (pj, m (r) = 

y/2JmWi,mr/R)/[RJm+iWj, m )], where j = 1,2,3,... 
and the argument f3j m is the jth zero of J m (x). More 

explicitly, we have ul ms (r) = T,j cl maj (f>j,m(r) and 
(r) — £\ d nmsj (pj tm+ 2(r) for the f radial wave func- 
tions, and ui, ns (r) = £\ ■■cimsj<t>i,™+\{r) and vi ms {r) = 

dtmsj ( f > j,rn,+i( r ) f° r the I ones, and they already satisfy 
the boundary conditions u° ms (R) = v° ms (R) — at the 
edge of the cylinder. 



Using 



the 



orthonormality 



condition 



/ rdr<pj vin (r)(j)ji >m (?") = Sjj> where Sjf is the Kro- 
necker delta, this procedure reduces the BdG equation given 
in Eq. (Q]i to a 4:j max x 4j max matrix eigenvalue problem, 
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A 33 ' 
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S 33 ' 
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—A JJ ' 












-A 33 ' 


- K 33 ' 


qj'j 
°m+l 




\ 


A 33 ' 





qj'j 


- K 33 ' 


) 



nmsj 1 \ 
'nrns-j' 
nmsj' 



\ nmsj' / 



< \ 

nrnsj \ 
nmsj 
nrnsj 

\ <f I 

\ nrnsj / 



(2) 



for each (m, s) sector, if we allow 1 < j < j max states. 
Here, K 3 a 3 ms = \0]J{2M a R?) + Pj(2M a ) - are 

the kinetic energy terms, S 33 = a jf 1 rdr4>j tTn (r)[d/dr + 
(m + l)/r](j>j' , m+ i(r) are the spin-orbit coupling terms 

leading to S 33 ' = aCy m f Q rdr<j>j, m {r)J m {^i tm+ \r/R) 

where C fm = \/2Pj, tm+1 /[R 2 J m+2 (0 jljm+1 )], and A 3 £ = 

J^rdrA{r)<j)j. m {r)<j)ji. m +i{r) are the pairing terms. The 
same procedure also reduces the order-parameter equation to 



A W = ^ \ C imsj d i.msj>h,ni+l{r)<t>3 , ,-m+2{r) 



nmsjj' 



X f{e nms ) + c n, ns jdn ms j'4'j,m( r ) < t > j' ;rn+l(r)f{ — Bums)}, 

(3) 



and the local-density equations to 

n t( r ) = TT-f E [ C nm«3 C lmsj'^i,m(r)^',m(r)/( £ nm a ) 



2ttL 

nmsjj' 

+ dl msj dl msj ,(l) jjm +i(r)(f>j', m+ 2( r )fi- E nms)}, (4) 

n ii r ) = E i C tnsj C imsj' < f > j,m+l(r)<Pj', m +l(r) 

nmsjj' 

x f(e nma ) +di msj df imsf ct> j , m +i(r)(t>j', m +i(r)f(-e rims )}. 

(5) 

We recall that the sums are only over the quasiparticle 
states with e nms > 0. Using the orthonormality condition, 
we also obtain the total number of a fermions as N a = 



F \(c° 



) 2 f{e nms ) + (d° msj ) 2 f(-e nms )}. We empha- 



size that these mean-field equations can be used for all values 
of a s and a at low T, but they provide only a qualitative de- 
scription of the system outside of the weak-coupling regime, 
i.e. in the BCS-BEC crossover. In this paper, we set the tem- 
perature to zero, and consider a strongly-interacting Fermi gas 
at unitarity, i.e. \a s \ —> oo, as a function of a. 



IV. NUMERICAL RESULTS 

In our numerical calculations, we set a large energy cut- 
off e c = lOep, and numerically solve the self-consistency 
Eqs. ©-(El) at T = 0. Here, e F = k 2 F /(2M) is a characteris- 
tic Fermi-energy scale where fejr is the Fermi momentum cor- 
responding to the bulk value of the total density of fermions, 
i.e. rif(r) + i%i(r) = fc|,/(37r 2 ) at the bulk. We also choose 
R = 25/kp as the radius and L = 10/kp as the length of the 
cylinder, and j ma x = 50 and |m| TOax = 75 as the maximum 
quantum numbers. Note that |s| TOaa; = L yj Me c /2 /it in order 
to be consistent with the energy cutoff. Since the presence of 
a single vortex line can not significantly effect the bulk param- 
eters, we first solve [i„ and Ao self-consistently for a vortex- 
free thermodynamic system, and then use these solutions as an 
input for our vortex-line calculation, where Ao corresponds to 
the bulk value of A(r). Here, we assume A(r) is real without 
losing generality. 



A. Order parameter and density of fermions 

In Fig. |2fa), we show typical order-parameter profiles A(r) 
for a = 0.5k F /M and a = k F /M when P = 0.5, i.e. 
it rapidly increases from zero around the vortex core, sat- 
urates to its bulk value Ao around kpr ~ 5 and then it 
rapidly decreases to zero near the edge of the cylinder ll27ll . 
Here, the population-imbalance parameter P = [nf(r) — 
n ^( r )]/[ n t( r ) + n ^( r )] i s defined at the bulk. We see that 
A(r) increases with increasing a, e.g. its bulk value increases 
from Q.50ep to 0.66ep, and that the effect of spin-orbit cou- 
pling is similar to the effect of increased interaction strength. 
This is due to the increased density of states with increasing 
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a, and it is consistent with the previous results on thermody- 
namic systems l3Hlql. 



(a) 



0.6 
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FIG. 2: (Color online) The order parameter A(r) (in units of ef) and 
density n a (r) (in units of k F /(3ir 2 )) profiles are shown as a function 
of radial distance r (in units of I/^f). Here, we set l/(fcfa s ) = 
and P = 0.5, and vary a. 



In Fig. |2jb), we show the corresponding density profiles 
n a (r) for the same parameters. When a = 0, it is well-known 
for the fermion superfluids that the density is depleted near the 
vortex core compared to its bulk value [27], and that the de- 
pletion increases with increasing interaction strength toward 
the molecular BEC side. This is because the energy separa- 
tion between the vortex core-bound states increases with in- 
creasing Ao which makes them less occupied. In fact, for a 
population-balanced Fermi gas, the density depletes fully and 
becomes zero at the vortex core in the molecular BEC limit, 
consistent with the theory of weakly interacting atomic Bose 
gases. However, for population-imbalanced Fermi gases, the 
vortex core may still be filled with excess fermions toward this 
limit (HI]. When a ^ 0, in Fig.[2|b), we again see that the ef- 
fect of spin-orbit coupling is similar to the effect of increased 
interaction strength, i.e. density depletion also increases for 
both a components with increasing a. To further understand 
the density depletions, next we analyze the spectrum of energy 
eigenvalues . 



B. Energy spectrum 

In Fig. [3] the spectra of energy eigenvalues e nms are shown 
as a function of planar angular momentum m for the s = 
sector. The spectrum rapidly becomes symmetric around 
rn = with increasing s, since the vortex core and edge states 
disappear when \s\ ~ 1. Here, we choose a and P such that 
Fig. [3 a) corresponds to a topologically-trivial gapped bulk SF 
phase, and Fig. Ob) corresponds to a topologically-nontrivial 
gapless bulk GSF phase (see the thermodynamic phase dia- 
gram given in Fig.Q]). First of all, we note that the excitation 
spectra shown in these figures have the necessary symmetry 
£nms = — £n,-{m+2).si which follows from the particle-hole 
symmetry of the Hamiltonian. In addition, a second branch 
of continuum spectra appears in both cases and on both posi- 
tive and negative energy regions when |e nms | > 1.2ef- This 
is similar to what happens in a thermodynamic system, for 
which the excitation spectrum has two quasiparticle and two 
quasihole branches when a ^ 
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FIG. 3: (Color online) The energy spectrum £„ ms (in units of ef) is 
shown as a function of planar angular momentum m when s = 0. 
Here, we set l/(kFa s ) = and a = 0.5&F /M, and vary P. In (a) 
we choose P = 0.2 corresponding to a topologically-trivial gapped 
bulk SF phase, and in (b) we choose P — 0.5 corresponding to a 
topologically-nontrivial gapless bulk GSF phase. The linear branch 
in (b) corresponds to the edge-bound states. 

When a = 0, the positive- and negative-energy spectra are 
connected by a single branch of discrete Andreev-like bound 
states 1271 |29fl. While the visible discreteness of the contin- 
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uum spectrum is a finite-size effect and the spectrum becomes 
continuous only in the thermodynamic limit (kpR — » oo), the 
discreteness of the bound states is insensitive to the system 
size since these states are strongly localized around the vortex 
core. The lowest-energy quasiparticle excitation requires a 
minigap of order Aq/(2£f ) <C An. When a^O, the energy 
spectrum of the SF phase shown in Fig. [3ja) is very similar to 
the usual population-balanced s-wave superfluids, for which 
the bulk energy spectrum is also gapped. There are only a few 
discrete core-bound states appearing within the bulk energy 
gap in Fig. |3ja), since the bulk order parameter A ~ Q.lep 
is quite large for P = 0.2 leading also to a large minigap. 
However, the bulk order parameter decreases to Ao ~ 0.5£f 
when P — 0.5, leading to a smaller minigap in comparison to 
P = 0.2 case, and hence a larger number of core-bound states 
as can be clearly seen in Fig. [3b). 

In contrast, we see a major difference in the energy spec- 
trum of the GSF phase as shown in Fig. [3jb). In addition to 
the branch of discrete core-bound states that is also present in 
the SF phase, there is a second branch of bound states which 
are strongly localized around the edge of the cylinder. These 
states result from Andreev scattering at the rigid walls of the 
cylinder, and their spectrum is linear in energy within the 
continuum gap [30]. We find that the lowest positive-energy 
and highest negative-energy bound states have m = — 1 and 
s = 0, and their energies are £n « 4.33 x 10~ 4 £f and Eq' ~ 
—4.33 x 10 _4 £f, respectively. This is not a coincidence since 
we know that the energy spectrum has e nms — — £ n ,-(m+2),s 
symmetry, and given that the spectrum is expected to have 
a two-fold degenerate zero-energy bound states, i.e. a pair 
of Majorana quasiparticles, in the thermodynamic limit, they 
must occur at rn = —1. However, hybridization between the 
core- and edge-bound states (see below) lifts this degeneracy 
in a finite system, and the zero-energy bound states split in 
energy as we find here. 

0.04 r- — , . , . — —, 




-0.04 1 1 1 1 1 1 

5 10 15 20 25 

k F r 

FIG. 4: (Color online) The radial wave functions w^ ms (r) and 
i)J ms (r) (in units of 1//cf) are shown as a function of radial dis- 
tance r (in units of 1 /fei?) for the lowest positive-energy bound state. 
For the parameters considered in Fig.[3lb), eo ~ 4.33 x 10 -4 £f and 
it occurs at m = —1 and s = 0. 

In Fig.|4] the radial wave functions u° ms (r) and v^ ms (r) 
of the lowest positive-energy bound state with energy £ are 



shown as a function of radial distance r, for the parameters 
considered in Fig. Ob). (Due to particle-hole symmetry, the 
radial wave functions for the highest negative-energy bound 
state with energy Eqi can simply be obtained by changing 
u nms( r ) v nms( r ) an ^ v i ce versa.) We note that while 
the wave functions have u^ ms (r) = — v^ ms (r) symmetry 
near the vortex core, they have u° ms (r) = v* ms (r) sym- 
metry near the edge. These are consistent with the symme- 
tries of Majorana quasiparticles IT26ll . for which the associ- 
ated quasiparticle operators are self-Hermitian, i.e. a Majo- 
rana quasiparticle is its own anti-quasiparticle. This is clearly 
seen from the Bogoliubov-Valatin quasiparticle creation oper- 
ator 7 t = f drJ2 a [K(r)ipi(r) + <(r)^ CT (r)] evaluated at 
in = — 1 and s = 0. Since the Majorana quasiparticles always 
come in pairs, they appear simultaneously but away from each 
other in real space. In our single vortex line, while one of 
them is mostly localized at the vortex core, the other one is 
mostly localized at the edge, with some degree of hybridiza- 
tion between them due to finite-size effects. The hybridiza- 
tion is clearly seen in the wave functions shown in Fig. [4] 
We note that due to this coupling between the Majorana core- 
and edge-bound states, their two-fold e nms — degeneracy 
is lifted, causing a small level splitting as discussed above. 
Increasing the separation between Majorana core- and edge- 
bound states, i.e. when kpR — > oo, weakens the hybridiza- 
tion such that both bound states eventually become degener- 
ate in energy with e nms = 0. When this happens, the core 
quasiparticle is well-localized around the vortex core with 
u nms( r ) = ~ v nms( r ) symmetry, and the edge quasiparticle 
is well-localized around the edge with u^ ms (r) — v^ ms (r) 
symmetry, without any hybridization between the two. 

So far, we have established a major difference between the 
energy spectra of SF and GSF phases, which is mainly due 
to the appearance of edge- and Majorana zero-energy bound 
states, and this difference leaves its signatures in various ob- 
servables as discussed next. 



C. Single-particle density of states 

For instance, the local single-particle density of a states 

D a (r, U ) = EJK(r)| 2 % - e n ) + K(r)| 2 <5(w + e n )], 
where S (x) is the delta function, as well as the integrated 
single-particle density of a states D a (uj) = J drD a (r, lu) 
provide direct evidences for the existence of edge- and Majo- 
rana zero-energy bound states as shown below. In particular, 
for a vortex line, and after using the orthonormality conditions 
for the Bessel functions, D a (uj) reduces to 

ArM - [«™^) 2(5 ( W ~ £ ™s) + {dlmsjf 5 ^ + £ nms}. 

nmsj 

(6) 

We use a small spectral broadening (O.OIef) to regularize the 
delta functions in our numerical calculations. 

In Fig. [5] the integrated density of a states D a (uj) are 
shown as a function of energy w, for the parameters consid- 
ered in Fig. [3] We see that while the density of states vanishes 
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(a) 



1400 



mitian conjugate. Since 3 a (r) circulates along the 
direction, i.e. J CT (r) = J a (r)0, we find 




FIG. 5: (Color online) The integrated density of a states D a (u) (in 
units of 1/ef) are shown as a function of energy uj (in units of ef), 
for the parameters considered in Fig. [3] 



around to = in the SF phase, due to the presence of a gap 
in the energy spectrum, it is finite in the GSF phase with very 
small peaks around uj — 0. These peaks are due to the pres- 
ence of discrete core and edge states within the continuum gap 
in energy, and they are most clearly seen in the majority (f) 
component. We also note that the appearance of a second con- 
tinuum branch in the excitation spectrum increases the density 
of states considerably when |e nms | > 1.2ef- Next, we ana- 
lyze the local probability-current density of a fermions, which 
also shows signatures for the edge- and Majorana zero-energy 
bound states. 



D. Probability-current density 

Similar to the usual a = treatment, the quantum me- 
chanical probability-current operator for a fermions can be 
identified from the continuity equation. While the presence 
of a spin-orbit coupling leads to additional terms in the total 
particle current operator, these terms do not contribute to the 
current since the expectation value (ipt(r)ip±(r)) = 0. There- 
fore, using the Bogoliubov-Valatin transformations, the local 
current density J CT (r) = [l/(2M,i)](^( r )V^(r) - H.c.) 
circulating around a single vortex line becomes 
J CT (r) " = [l/(2M CT i)]£ n K*(r)V<(r)/( £ „) + 

<*(r)V<(r)/(-£ n ) - H.c], where H.c. is the Her- 



J t (r) 



2lT ~M^ Yj { m E C L SJ ^,mW] 2 /(enm S ) 

' nms j 

(m + 2)[^4 msj J , m+2 (r)] 2 /(-£ nms )}, (7) 

3 

)] 2 /( 



2TrM±r 



(m + l)[Y,di msj ^ m+1 (r)} 2 f(-e nms )}, (8) 



for the strengths of the particle-current densities. 
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FIG. 6: (Color online) The probability current density of a fermions 
(in units of M/kp) are shown as a function of radial distance r (in 
units of 1/fcF), for the parameters considered in Fig. [2] 

In Fig.|6l the probability-current density of a fermions are 
shown as a function of radial distance r, for the parameters 
considered in Fig. [2] When a = 0, the core-bound states 
have negative (diamagnetic) and the continuum states have 
positive (paramagnetic) contribution to J CT (r). This leads to 
a nonmonotonic J<r(r) which first increases as oc r and then 
decreases as oc 1/r 12711 . The latter behavior is due to the sat- 
uration of the superfluid density for long distances away from 
the vortex core. Therefore, a maximum peak current occurs 
at some distance away from the vortex core. When a ^ 0, 
the major difference is at the edge. The j" and J, currents flow 
in opposite directions, and their magnitude increases with in- 
creasing a. Such a counterflow of mass currents occurs even 
for the SF phase (not shown). Since the sums in Eqs. (0 
and ([8]) are over states with e nms > 0, counterflowing edge 
currents result from the asymmetry of the energy spectrum 
around m = 0, due to the presence of edge states. 



V. CONCLUSIONS 

To conclude, we considered a Rashba-type spin-orbit cou- 
pled Fermi gas, and used the BdG formalism to analyze a sin- 
gle vortex line along a finite cylinder with a periodic boundary 
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condition. When the populations of the pseudo-spin compo- 
nents are sufficiently imbalanced, depending on the strength 
of the spin-orbit coupling, we showed that core- and edge- 
bound states as well as Majorana zero-energy quasiparticles 
appear in the energy spectrum. These states leave signatures 
in the density of single-particle states and particle-current den- 
sity, and particularly, we found that the pseudo-spin compo- 
nents counterfiow near the edge of the cylinder, the strength 
of which increases with increasing spin-orbit coupling. 

While preparing this work, we became aware of a closely 
related work lf3ltl . where the vortex core and edge states are 
analyzed for a trapped two-dimensional Fermi gas. For the 
most parts, our work is consistent with their findings. How- 
ever, in contrast to our finite-cylinder setup, where the system 
is either in an SF or a GSF phase, depending on the parameter 



regime, both SF and GSF phases may also coexist in a trap 
in different regions. The possibility of such a phase coexis- 
tence again leads to Andreev scattering at the SF-GSF phase 
boundary, giving rise to an additional branch of edge-bound 
states. 
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